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A fUQDIPIED DODGE AlGORITIfM FOR THE PARABOLIZED NAVIER-STOKEvS 
EQUATIONS AND COMPRESSIBLE DUCT FLOWS 

By 

C. H. Cooke* 

SUMMARY 

A revised version of Dodge’s split-velocity method for numerical cal- 
culation of compressible duct flow has been developed. The revision incorp- 
orates balancing of mass flow rates on each marching step In order to main- 
tain front-to-back continuity during the calculation. Tlte (checkerboard) 
zebra algorithm is applied to solution of the three-dimensional continuity 
equation in conservative form. A second-order A-atable linear multistep 
method is employed in effecting a marching solution of the parabolized 
momentum equations. A checkerboard iteration is used to solve the resulting 
implicit nonlinear systems of finite-difference equations which govern step- 
wise transition. Qualitive agreement with analytical predictions and exper- 
imental results has been obtained for some flows with well-known solutions. 


*Professor» Department of Mathematical Sciences, Old Dominion University, 
Norfolk, Virginia 23508. 
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INTRODUCTION 

It has been said that the full Navler-Stokes equations represent the 
ultimate mathematical model upon which to base numerical algorithms for 
predicting flows of practical significance! However » even with the advent 
of the so-called vector computers with vast virtual memory and quadrupled 
processing speeds » extant numerical and computational difficulties are 
sufficient to merit a search for simpler mathematical models and less com- 
plicated numerical methods which can still provide useful solutions to prob- 
lems of Interest* Thus, considerable analysis and numerical experiment has 
been devoted to the exploitation of parabolized marching methods for flow 
prediction. References I to 7 ^represent a perhaps typical but by no means 
exhaustive sampling of the available literature on this subject. 

The parabolized marching methods are somewhat more general In applica- 
tion than the classical boundary-layer approach, since transverse pressure 
gradients are not disregarded and, In some cases, upstream Influences can be 
transmitted through the pressure field. However, the basic assumption that 
stroamwise viscoua diffusion can he neglected restricts application £0 flows 
with a primary flow direction, limited upstream influence, and which may 
exhibit, at worst, crossplane recirculation* Unfortunately, in subsonic and 
transonic wind-tunnel flows, the elliptic upstream Influence can be a slg*^ 
niflcant factor in the flow dynamics; hence, Interest arises in simpler 
mathematical models which permit this interaction, A case In point has been 
the development of Dodge's velocity splitting method, which allows global 
propagation of influence through the pressure field and which has met with 
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successes In both unconfined compressible and confined compressible flows 
(refs. 7*10). However* the method is not yet fully proven. 

In this paper we shall be concerned with the application of a compres* 
sible formulation of Dodge's split velocity technique (ref. 9) to the calcu- 
lation of developing flow in a square duct. The original method has been 
revised to affect constant mass flow rate on each transverse plane while 
marching down the channel. Parabolized momentum equations are employed. 
However* a fully elliptic pressure field is allowed by the Iterative manner 
in which the solution of the continuity equation is coupled into the calcu- 
lation procedure. Application of the presently developed computer algorithm 
is restricted to subsonic flow. It could readily be altered to allow tran- 
sonic calculations through modification or replacement of the algorithm used 
to solve the conservative continuity equation. 

Computational simplicity as well as numerical stability is achieved in 
marching the momentum equations with an A~stable (ref. 11) implicit linear 
multistep method* the equations of which are iteratively solved at each step 
by employing checkerboard successive overrelaxation. Wlille this solution 
procedure may be considered expensive, the presence of quadratic as well as 
higher order nonlinearities in the parabolized momentum equations requires 
that some iteration be employed to Improve accuracy. As an extra benefit, 
the wide-ranging stability of the resulting marching equations appears well 
worth the cost. 

Finally, the peak efficiency of the methods developed is undoubtably 
best realized on the computer system for which it is has been designed* 
namely* the Cyber 203. For puch machines, a numerical algorithm must 
effectively exploit the array-processing capabilities; otherwise, methods 
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which ate not highly vectotitablo ralauae the available computing potential 
and can result in quite ordinary processing speeds* The explicit nature of 
the checkerboard algorithm yields a highly vectorlrable method ideally suit- 
ed for tlw array processor. 

In certain parabolized marching schemes for confined flow (ref. 1), it 
has been the practice to decouple streamwlse and transverse pressure gradi- 
ents. Some argue (ref. 12) that this is necessary in order to obtain mean- 
ingful physical solutions with parabolized equations. While results are 
still Inconclusive, computational experience gained in the current research 
appears to support this belief. .Weak, but not total, uncoupling of the 
streamwlse pressure gradient has appeared necessary, altho«»sh this may stem 
from the manner in which local ^/ontinulty of mass flow Is enforced. 

As noted by Pa tanker and Spalding (ref. 6), the description of a numer- 
ical procedure for solving the Navler-Stokes equations can have two aims, 
which ate seldom possible or desirable to accomplish simultaneously. The 
first aim is to convey to the reader an understanding of the major prin- 
ciples in sufficient detail that someone with a background in the area could 
Improvise the remainder for himself. The second aim is to present the par- 
ticular equations and all approximations employed to a degree that the com- 
putational experiment could be Identically reproduced. However, the second 
mode requires such proliferation of detail that smooth reading is Impeded 
and understanding is inhibited. Tlierefore, the first aim has been chosen 
for the present paper, and this will be attempted in the following sec- 
tions* 
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LIST OF SYMBOLS 


specific h^ats 
static pressure 

I 

density 

tli/ree" dimensional velocity vector 
■vyiscosity 

Reynolds numberi ^o^o^ 

C/C 

p' V 

temperature 

reservoir values for tet^erature, pressure, density, 
speed of sound, and viscosity 

acalar potential 

relaxation of parameter 

channel half~width 

Mach nunber 

ratio of channel streaniMise velocity component to 
entrance value 


PARABOLIZED GOVERNING EQUATIONS 


The nondJLtnenslonal Navier-Stokes equations for compressible steady flow 
with which we shall be concerned are presented below* 

Continuity; 

7»pw“0 (1) 


Momentum: 


p(w • V)w ■ ~VP 


Energy; 

1 “ 2 

T » T_ - w 
° 2 


7x ~7xW 
Re 


+ 7 


! l . JL. 

3 Re 


w 


( 2 ) 




Hare> for flow in ducts with nonconduGfclng walls j the usual energy equation 
has been replaced by the algebraic constant total temperature relation [eq. 
(3)]. The constitutive relations are 
Y ~ 1 


P 


PT 


(4) 


and the viscosity approximation 

p - (y - i)t. 


(5) 


For subsonic flow the governing equations are elliptic. However, a 
common approximation used to parabolize these equations (refs. 1,2) is ob- 
tained by neglecting streamwlse diffusion terms in equation (2). With the 
exception of the entry region, the approximation is considered valid for 
flow in channels whose lengths are large compared to half-width (ref. 2). 
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Perhaps it should be remarked that when Dodge's method is applied in ob~ 
taining numerical solutions of these equations the approximation is only a 
partial parabolization, since the pressure field is obtained from an 
elliptic boundary value problem. This, of course, allows global propagation 
of disturbances, through the pressure field and the iteration process. 

a 





DODGE’S HETHOG 


Introduction 


In Dodge’s method, the total velocity vector w is arbitrarily de- 
composed as a sum of rotational and irrotational parts; 

w »• d* u (6) 

where ^ is a scalar potential, Pressure is hypothesized to depend solely 
upon the irrotational velocity according to the Isentropic relationship 

However, density is decomposed as the am of a viscous contribution 
and an isenttopic contribution p*: 


P - Py ^ p’ 


ri 


where 


P* ■ P, 


1 - 


oi ai 


l/Y-l 


Substituting equations (6), (7), and (9) in equations (1) and (2) leads to 
the so-called split equations ^ 

V • p7^ " - V • pu (10) 


u „.r 1 _ nM u 


P(w • 7)w - p*(V4» • V)Vt + Vx rg- 7xw - V V • w U 0 (U) 


J V" 7 

Equations (10) and (11) are to be iteratively solved: equation (10) with a 
3-D relaxation method for elliptic equations, following which the para- 
bolized version of equation (11) is marched downstream by employing a check- 
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erboard Iteration to aolve an Implicit gystem oi equations at each step. A 
synopsis of the iteration procedure la now presented* 

Overview of Iteration Procedure 

(1) Determine a suitable initial pressure distribution by estimating 
a global ^ distribution* In this investigation! pressure on the 
first pass is assumed U be a function of only streamwlse displacement^ 
and a mas8**balancing operation estaMishes the initial pressure field* 

(2) Employing the current pressure field, march a parabolized version of 
equation (11) down the duct, simultaneously storing the right'“hand side 
of equation (10). (See also eq. (17)]. 

(3) Solve equation (10) (or eq. (17)] to obtain an updated pressure field. 

(4) Repeat the computational pass consisting of steps (2) and (3) until 
sufficient passes and a converged pressure field are obtained* 

Dodge's Method Revised 

Dodge (ref. 9) reports problems arising from adjustment of front^to- 
back continuity requlremextts with an iteration which is similar to that 
previously outlined. It is expected that this slow convergence stems from 
Incomplete satisfaction of the continuity equation which could, for example, 
be solved after the momentum march terminates in some form such as 

V • p - - 7 • (pu) , (12) 

'^n ^n n-1 

This is in contrast to the usual parabolized marching methods, for which’ 
both mass and velocity variables are updated at each marching step. 

Physically, in order to maintain continuity in a channel flow, the mass 


« * //pfll + uVydz - 


normal 


(13) 


must remain constant at each transverse plane* However, in Dodge’s (unre- 
vised) method this provision is only weakly incorporated through equation 
(10), which is solved globally upon termination of a marching pass* Thus, 
poor satisfaction of mass balancing during the momentum marching process is 
only to be expected, as numerical experimentation indicates* 

Consequently, we have chosen to revise the Dodge technique in a manner 
which alleviates this difficulty* This was at first attempted by employing 

^ - I"" s(K)dK + ?(x,y,z) (14) 

to write equation (12) in the form 
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V . (pv,)^ - - V • (pu)„_, - 


(15) 


The function g is determined by iterating the numerical counterpart 
of the parabolized equation (11) at each fixed marching step until numerical 
balance of mass flow rate is achieved* This is accomplished through gradual 
changes in streamwise velocity, pressure, and density effected by the equa- 
tion 


k+1 




(16) 


normal' 

with a a relaxation parameter. Aside from the benefit of an Instanteous 
balance in mass flow rate, another merit of this device is that fewer global 
Iterations are required in the relaxation solution of equation (15), as it 
is now more nearly satisfied at the outset* 
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However, this approach was found defective, In theory as well as in 
fact* The solving of equation (15) in the form indicated yields nonphysical 
results, as it provides a quasi-full potential transonic flow equation whose 
elliptic-hyperbolic transition point can differ markedly from Mach I* This 
difficulty can be largely alleviated, although not totally circumvented, by 
replacing equation (15) with the equation 

Sf • - V * (pu)^^j (17) 

whose point of transition more closely approximates the physics of the flow. 
Pressure gradients in Dodge’s unrevlsed method would be computed on 


pass n from the equation 

r 






/ 1 o\ 
\ *a/ 


In the revised version, pressure gradients are allowed to develop 

« 

during the raass-balancitig Iteration according to the equations 


9x 


• Pk.n'S Sx + S ^xy *xz> 


<19) 


- - p! (gS” + 9" + <>" +" ) 

8y *^k,n'-® ^xy ^y ^yy ^yz-' 


3P 

9z 


n,k 


^xz ^y ^yz ^z ^zz-' 


( 20 ) 

( 21 ) 


The quantity g is determined through equation (16), and by 

second-order backward differencing. This precedure represents a weak de- 
coupling of the streamw#'?c pressure gradient, since the g terms are the 
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dominant cont;;ibution8, and since these contributions are determined from 
local plane- to-plane continuity considerations, somewhat independently of 
the output from the global continuity equation [eq. (17)1 on the previous 
pass. 
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NUMKRICAt mVtBU 


Thck algorithm deemed moat efficient for munerically solving equation 
(17) on the atray-proceaaing computer ia the tebra olgoritbn of South et al. 
(ref* 13). Tltia 3“t) relaxation technique la in some reapectB aimllar to the 
hopacotch method of Gourlay (tel. lA)* In equation (17) central differencea 
are applied to all derivative terma. Variables in plane i are updated in 
checkerboard faahiont plane by plane in a downstream aweep) using already 
updated values at plane 1 ** I and old iteration values in plane i + U 
tterativa repetition of downstream aweepa la used to converge the field, 
with a relaxation parameter employed to speed convergence. 

A second-order accurate, implicit linear multistep method is used on 
equation (11) to march in the stepwise direction* Ihe implicit equations 
are Itoratlvaly aolved using a checkerboard successive overrelaxation 
scheme, with maaa hnlnncing built in as previously described. Streamwise 
darivatlves ate backward differenced second-order accurate, while deriv- 
atives in the (transverse) cross'^plane are approximated second-order using 
central dlfftrences. A prediction of form 

^l ** ^ 1-2 

la used to initially estimate a velocity variable in plane 1. Hie checker- 
board method is then employed on the differenced counterpart of equation 
(11) to update variables in plane i in two cycles, with values updated on 
cycle 1 fad into the succeeding cycle. Tt\ls two-cycle update process is 
iterated, employing equations (16) and (19) to (21) to alter the flow speed 
and pressure gradients until a balance in mass flow is achieved. 
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DEVELOPING FLOW IN A STRAIGHT DUCT 


Tha revl^iod method of Dodge haa been employed to develop a flnite-dlf~ 
ferenoe numerical model for three-dlraenaional vlacoua flows in confined 
regions. For bound ary-layer resolution, the capability to allow individual 
coordinate stretching in each coordinate direction has been Incorporated. 

The method so developed has been programmed using the SLI vector language 
for the Cyber 203 array processor, and appears debugged. Tlie 32-bit half- 
word option of SLI has been employed in programming the aebra relaxation 
algocittwn for solving equation (17), while 64-bit full-word arithmetic is 
used in programndug the checkerboard marching algoritlwi. live program has 
boon tested by application to tho problem of computing the steady developing 
flow in a straight duct (see fig. 1). Boundary conditions for tha problem 
are now given. 

Boundary Conditions 

■ Inflowl T - T - -i u7 

° 2 i 

WM mm 

specified velocity profiles, W^ •« H(y,a), 

“ R(y»55)* •' constant 

'l'j^(O.ypa) • 8(0) 8 q 

Duct walls; velocity no slip, 4 “ 0, T • T , p •" p 

11 W W 

Outflow; p^ extrapolated, (Ji g^ extrapolated 

Artificial barriers; Tl\e computational domain is taken to be one 
quarter of the total duct cross section, and symmetry conditions 
are applied at the two resulting (nonwall) artificial barriers. 

Hare the normal velocity component vanishes together with normal 
derivatives of «(i and the other velocity variables. The variables 
P and T, of course, depend on and velocity at these boundaries. 
However, for constant total temperature, vanishing normal derivative 
in T, p is the natural boundary condition. 
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COMPUTATIONAL RESULTS 


'iv 

Introduction 

An asBQSsfflont of approximation error Inheront in numarlcal solutions of 
parabolized flow models computed with the revised Dodge's method has been 
undertaken* Ttiis has been accomplished for certain flows wltose solutions 
have been either analytically or experimentally determined and which are 
available for comparison with numerical results. 


Two-Dimensional Channel Plow, Re •« 100 
Numerical solutions for two-dimensional, low Mach number (M > 0.05) 
flow in a straight channel have been computed, employing a uniform mesh of 
17 X 17 X iQO nodes. The development of the normalized centerline velocity 
component Vj^/U characterizing the two-dimensional channel imd the corres- 
ponding results of Goldstein and Kreid (ref. 15) for three-dimensional flow 

* 

in a square duct are compared in figure 1. For the fully developed case 
Schllcting (ref. 16) gives a limiting value Vj^/U ■ 1.52, Tlie corres- 
ponding maximum value from the computational results of figure I is - 
1.A25, This represents a relative error of six percent or less, depending 
upon how fully developed the numerically calculated flow is considered to 
be. One should perhaps also bear in mind that Schlicting’s results come 
from matched asymptotic expansions, which are in themselves approximate 
methods. Figure 2 exhibits typical spanwise deviation from a parabolic 
profile characteristic of at nimierical full development. The maximum 
error is again on the order of five or six percent. 

15 
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Three Dimensional. Duct Flow, Re « 100 

i, 

A 32 X 32 X iOO rneiuh has bean employed to compute low ^tach number* 
threa-dimanf^lonal flow in a square duct. Straamwise variation in centerline 
velocity component and pressure F are shown Iq figures 3 and 4. A 
comparison between numerical calculations for V^/Vi and the corresponding 
experimental results of Goldstein and Kreld (ref. 15) is shown in figure 5. 
Comparative fully developed values for Vj^/O at centerline for the two 
cases are 1.93 and 2.10. Including grid stretdilng in the calculation im- 
proved the former value to 1.9725. Tl\i8 again represents a relative error 
in the vicinity of six percent. 

Figures 6 to 8 show numerical performance Indicators which trace the 
convergence history of a typical calculation. Here two kinds of outflow 
boundary conditions for solution of the elliptic velocity potential equation 
have been tested. For the calculation represented by figure 8* a Dlrichlet 
condition was applied. Although a seemingly more rapid* nonosclllatory 
convergence history is observed* distortion in straamwise velocity and pres- 
sure profiles near the channel exit accompanied its application. A Neumann 
boundary condition appeared to eliminate these nonphysical distortions* but 
at the misfortune of a more sluggish* oscillatory, convergence history (see 
figs. 6-7). The Neumann condition was chosen for the calculations hereto- 
fore discussed. 
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Figures 9 and 10 show characteristics of the exit crossflow and bound 
ary-layer profiles for the streamwlse velocity component. In calculating 
the fully developed crossflow of figure 10, a channel length of * 12 
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units with a coavsa mash having 144 downstraam nodes was employed. The sink 
itregularity near channel centerline in figure 9 was observed to move closer 
to channel center as the length increasedi and appeared to have coalesced 
with lb in the calculation of figure 10. The tendency of the core flow to 
be directed towards channel centerline (see figs. 9 and 10) is in qualita- 
tive agreement with the theoretical predictions of Rubin (ref. 4) concerning 
the asymptotic behavior of duct crossflow. 

CPU Time 

A computer code for the revised version of Dodge's method employs the 
Shl~language for the Cyber 203 version of the CDC-STAR computer. With a 
17 X 17 X 196 mesh, measured time for program execution was found to average 
2.3 X 10”*' computer resource units (CRU's) per node per model equation per 
pass. (This aasumss only four model aquationa. since the energy equation 
was replaced with an algebraic relationship.) A pass consisted of marching 
the momentum equations once down the duct, iteratively solving an implicit 
system for velocity propagation at each streamwlse station using a checker- 
board scheme, and Chen finding the zebra solution of the three-dimensional 
velocity potential equation* When measuring performance in CPU seconds per 
node per model equation per pass, the corresponding figure becomes 4.14 x 
lO”**. 
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SUMMARY AND CONCLUSIONS 


A revised version of Dodge’s velocity-split: method for numerical solu- 
tion of compressible confined flow has been developed* Numerical results 
from test calculations for low Reynolds number flow appear encouraging* In 
particular, qualitative but not wholly satisfactory quantitative agreement 
between present calculations and analytical predictions of Schllctlng (ref* 
16) and Rubin (ref. 4), together with experimental measurements of Goldstein 
and Kreid (ref* 15), has been achieved* 

However, the method la by no means fully understood or exhaustively and 
conclusively tested* A curious feature of the present approach Is the need 
for weak decoupling of the strearawlse pressure gradient, found necessary In 
order to achieve a convergent numerical solution* However, Spalding (ref* 
12) alleges that to obtain meaningful solutions using parabolized equations 
such a full decoupling is necessary* Further, Briley (ref* 1) reports suc- 
cessful and quantitatively accurate calculations obtained with an algoritlmi 
which incorporates this practice* Be that as It may, in this study conver- 
gence problems, first manifested by irregularities In the entrance region of 
the duct, were observed, but disappeared upon weakly decoupling the stream- 
wise pressure gradient. Even so, the numerical calculations appeared overly 
sensitive to the treatment of the Inflow, a characteristic of parabolized 
duct flows perhaps hinted at by Briley’s suggestion of shutting off the 
crossflow in the first few steps of the march (ref* 1). Rubin (ref. 4) 
alleges that the full Navier-Stokos equations are required to properly model 
entry region flow. Hence, the failure to properly smear errors in this 
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region may have been the root problem, assuming such errors can be widely 
propagated, either upstream or down. In solving the ellljitlc velocity poten- 
tial equation* 

In addlton to cases previously recorded, a convergent numerical solu- 
tion for Re - 1,000 flow in a duct of length ^ " 6 has been achieved. 
However, for satisfactory calculations at high Reynolds numbers, or for 
longer ducts, a more powerful adaptive grid capability than is to be afford- 
ed by simple stretchings in individual coordinates appears necessary. In 
duct flow the driving mechanism Is streamwise boundary-layer growth and 
concomitant flow acceleration in the core, with the induced crossflow 
strongest away from the walls. Hence, as the flow develops near-wall grad- 
ients decrease, or spread In extent, while core gradients grow, lor high 

Reynolds numbers the result is that a grid system with a simple stretching 
* 

mechanism set to capture entry region boundary layers is likely to become 
insufficient for resolving the more complex global pacterns which emerge as 
the flow develops. Since boundary-layer thickness Is thought to grow pro- 
portional to (x/Re)^'^^, It appears that the ideal grid should adaptively 
relax near-wall clustering in some such fashion. 

The revision of Dodge's method reported herein as regards mass balanc- 
ing is new, although classical In its physical origins and used previously 
with many other computational schemes. For flow In a straight duct It 
appears to be highly useful. However, thv'i original intent of this invest- 
igation of Dodge's method was to discover its possible utility as a tool for 
solving the slotted wind-tunnel problem, where it was felt that the upstream 
Influence permitted by the elliptic velocity potential equation would be 
highly desirable. This could well be so. Were the method found workable 
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without resort to the mas a- balancing stratagem. However, it is not telt 
that the present forroulation with added on mass-balancing can be useful. 

This is because of the possiblity of large pressure variations from tunnel 
to plenum, which would make the simultaneous balancing of mass on a plane 
extending over two essentially disparate channels with limited communication 
through the eslots unfeasible, as balancing over the entire region would have 
to be accomplished through tuning (the streamwise pressure gradient. 

Possibly the roost successful innovation of this investigati^m is the 

use of the checkerboard iteration as a tool for solving at each marching 

plane the implicit f inite-differenced momentum equations. It should be 

noted that the checkerboard algorithm used is different from either the 

hopscotch method used by Rudy (ref. 17) or the recent hopscotch innovation 

of Greenberg (ref. 18). In both cases, iteration is not practiced in 

solving an implicit system; hence, the well-known inconsistency of xihe 

hopscotch method becomes significant, at least during the transient calcu- 
li 

i! 

lation. In the present case no transient exists, since the steady equations I 

are being solved. The inconsistency of the usual hopscotch method is thus 
hard to bear, as it could lead to problems in solution accuracy. 
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Figure S. Exporlrflental veraua computed development for atreamwiae 
velocity component - throe-dimenaional duct flow. 
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Figure 7. Iteration, history, streamwise increment in velocity potential - 
Neumann condition at outflow. 
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Figure 8. Iteration history, corner and centerline pressure at outflow 
Dirichlet condition on velocity potential. 
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Figure 9. Crossflow and boundary-layer profile for streamwise velocity 
component at channel exit, Re » 100, ^ 3.0. 
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